Данная утилита реализует решение прямой задачи на геопотенциальные поля (поля силы тяжести и магнитные поля) - расчёт двумерного поля на регулярной 2D сетке на плоскости, исходя из плотностной или магнитной сетевой (сеточной) 3D модели среды. Задача - корректная и решается однозначно.

Обратите внимание, что, по сравнению с обычно принятыми формулировками, в ГИС INTEGRO ось z направлена вверх, но в пользовательском интерфейсе утилиты глубина задаётся в виде положительных значений, а система инвертирует их знак.

Математическая формулировка

Метод основан на численном интегрировании функции физического параметра, заданной на 3D сетке.

Формулировка моделей

Модель среды задаётся регулярной 3D сеткой ${\displaystyle \mathrm{{\rm Y}}=(x,y,z)\in {\mathbb{ℝ}}^3}$. Каждой точке сетки ставится в соответствие ячейка - прямоугольная призма, покрывая всю область моделирования без зазоров и перекрытий. Ячейки одинаковые (конгруэнтные), пересекаются не более чем своими гранями (с нулевым объёмом), и заполняют всё пространство между точками. Таким образом, модель является конечноэлементной с элементами - прямоугольными призмами.

Каждой ячейке ставится в соответствие одно значение физического параметра (плотности, намагниченности, восприимчивости...), и все точки одной ячейки имеют одно и то же значение физического параметра. Функция ${\displaystyle g(x,y,z),(x,y,z)\in {\mathbb{ℝ}}^3\wedge z<z_0}$ определяет значения физического параметра среды в области, охватываемой ячейками.

Модель поля задаётся как функция ${\displaystyle u(x_0,y_0)}$ от точек поля ${\displaystyle (x_0,y_0)\in \mathrm{\Omega }\subset {\mathbb{ℝ}}^2}$, ${\displaystyle \mathrm{\Omega }}$ - сетка поля. При этом в 3D пространстве все точки находятся на заданной высоте ${\displaystyle z_0}$, которая для простоты изложения ниже принимается за 0 без потери общности.

Предполагается, что сетки ${\displaystyle \mathrm{\Omega }}$ и ${\displaystyle \mathrm{{\rm Y}}}$ соответствуют друг другу - шаг, начало отсчёта в проекции ${\displaystyle Oxy}$ и азимут совпадают. По-другому: для каждой точки ${\displaystyle (x,y,z)\in \mathrm{{\rm Y}}}$ есть точка ${\displaystyle (x,y)\in \mathrm{\Omega }}$, и наоборот, ${\displaystyle (x,y)\in \mathrm{\Omega }⟹\mathrm{\exists }z:(x,y,z)\in \mathrm{{\rm Y}}}$.

Взаимосвязь между моделями

Известная взаимосвязь между распределением физического параметра в полупространстве ${\displaystyle V}$ под плоскостью ${\displaystyle z_0}$ и полем, заданным на горизонтальной плоскости:

${\displaystyle u(x_0,y_0)=\left[Ag\right](x_0,y_0)=\underset{V}{\iiint }g(x,y,z)k(x-x_0,y-y_0,z)dxdydz}$, (1)

где ${\displaystyle A=A(g)}$ - (линейный) оператор прямой задачи, а ${\displaystyle k(x,y,z)}$ - поле от точечной массы или магнитного диполя.

Гравитационная задача

Для поля силы тяжести определяется, что поле - первая вертикальная производная гравитационного потенциала: ${\displaystyle k(x,y,z)=\gamma \frac{-z}{r^3}}$, где ${\displaystyle r=\sqrt{x^2+y^2+z^2}}$; ${\displaystyle \gamma }$ - коэффициент, связанный с гравитационной постоянной на основе CODATA (на момент написания документации - CODATA2014). Гравитационная постоянная ${\displaystyle {\gamma }_{SI}=6.67408\times 10^{-11}{m}^{3}k{g}^{-1}{c}^{-1}}$. Учитывая единицы плотности ${\displaystyle g/c{m}^3}$, расстояния ${\displaystyle m}$, объёма ${\displaystyle m^3}$ и силы ${\displaystyle mGal}$, коэффициент ${\displaystyle \gamma =6.67408\times 10^{-3}mGal\times c{m}^3\times g^{-1}\times m^{-1}}$ или ${\displaystyle \gamma =6.67408\times 10^{-3}mGal\times m^2\times tonn{e}^{-1}}$.

Магнитная задача

Подразумевается, что данные измерений магнитного поля - скалярные, и измеряется проекция магнитной индукции на направление внешнего намагничивающего поля Земли. Тогда ${\displaystyle k=k(x,y,z)}$ состоит из произведения:

1. скалярной компоненты магнитной индукции, создаваемой веществом с единичной намагниченностью (магнитного момента единицы объёма) ${\displaystyle \overrightarrow{M}}$ в точке ${\displaystyle (x,y,z)}$;
2. ${\displaystyle C}$ - коэффициента для единиц размерности, ${\displaystyle C=0.1}$ или ${\displaystyle C=100}$ (в зависимости от единиц намагниченности).

Единицы измерения - ${\displaystyle \text{нТл}}$ для магнитной индукции, и ${\displaystyle \text{мА}\cdot {\text{м}}^{-1}}$ или ${\displaystyle \text{А}\cdot {\text{м}}^{-1}}$ - для намагниченности (в старых версиях ГИС ИНГТЕГРО поддерживалась единица намагниченности только ${\displaystyle \text{мА}\cdot {\text{м}}^{-1}}$).

В случае использования магнитной восприимчивости ${\displaystyle \kappa }$ она пересчитывается в намагниченность, используя усреднённый вектор индукции внешнего магнитного поля ${\displaystyle {\overrightarrow{B}}_{ext}}$: ${\displaystyle \overrightarrow{M}\approx \kappa {\displaystyle \frac{{\overrightarrow{B}}_{ext}}{{\mu }_0}}\approx \kappa {\displaystyle \frac{{\overrightarrow{B}}_{ext}}{4\pi 10^{-7}}}}$. Значения магнитной восприимчивости задаются в единицах ${\displaystyle 10^{-5}}$ СИ. При этом в формуле для ${\displaystyle k}$ используется константа ${\displaystyle C=100}$. Поддерживается вычисление ${\displaystyle {\overrightarrow{B}}_{ext}}$ из табличных данных (IGRF).

Дискретизация

Так как значение физического параметра в рамках элемента - прямоугольной призмы - не меняется, интеграл в (1) можно заменить на сумму интегралов по элементам:

${\displaystyle u(x_0,y_0)=\underset{V}{\iiint }g(x,y,z)k(x-x_0,y-y_0,z)dxdydz=\sum _{v\in {𝒱}}\underset{v}{\iiint }{g}_{v}k(x-x_0,y-y_0,z)dxdydz=\sum _{v\in {𝒱}}{g}_v\underset{v}{\iiint }k(x-x_0,y-y_0,z)dxdydz}$, (3)

где ${\displaystyle g_v}$ - значение физического параметра в элементе ${\displaystyle v}$, ${\displaystyle \mathrm{\forall }(x,y,z)\in v:g(x,y,z)=g_v}$.

Формулировка через дискретные преобразования Фурье

Если принять ${\displaystyle \mathrm{\Delta }X,\mathrm{\Delta }Y}$ - размеры ячеек по горизонтали, то можно заметить, что в (3) интегралы ${\displaystyle \underset{v}{\iiint }k(x-x_0,y-y_0,z)dxdydz=k_v}$ для ячеек на одной глубине отличаются только сдвигом. Это поможет ускорить вычисления с помощью быстрых преобразований Фурье.

Перепишем (3) через двойную сумму: 1) по ячейкам внутри одного горизонтального слоя ${\displaystyle {{V}}_{{𝒾}}}$ толщиной в один элемент и 2) по горизонтальным слоям:

${\displaystyle \sum _{v\in {𝒱}}{g}_v\underset{v}{\iiint }k(x-x_0,y-y_0,z)dxdydz={\displaystyle \sum _{i=1}^{n_z}}\sum _{v\in {{V}}_{{𝒾}}}{g}_v\underset{v}{\iiint }k(x-x_0,y-y_0,z)dxdydz}$. (4)

Рассмотрим тройной интеграл. Он зависит только от объёма интегрирования ${\displaystyle v}$, обозначим его как ${\displaystyle k_v(x0,y0)=\underset{v}{\iiint }k(x-x_0,y-y_0,z)dxdydz}$. Заметим, что ${\displaystyle k_v}$ для различных ячеек ${\displaystyle v}$ отличаются только сдвигом по координатам ${\displaystyle x_0}$ и ${\displaystyle y_0}$.

Теперь рассмотрим внутреннюю сумму. Она является расчётом линейной задачи от одного слоя элементов, соответствующих точкам сетки ${\displaystyle \mathrm{{\rm Y}}}$ с ${\displaystyle i}$-ым значением координаты ${\displaystyle z}$. Её можно представить как линейный оператор ${\displaystyle u_i=A_i{\hat{g}}_i}$ прямой задачи для горизонтального слоя толщиной в одну ячейку; аргументом оператора является вектор значений физического параметра ${\displaystyle {\hat{g}}_i\in {\mathbb{ℝ}}^{N({{V}}_{{𝒾}}{)}}}$. Так как отдельные ${\displaystyle k_v}$ отличаются только сдвигом, то в (4) внутренняя сумма - это непериодическая дискретная свёртка, а ${\displaystyle A_i}$ - матрица Тёплица (в 2D). Аппроксимируем их периодической дискретной свёрткой и циркулятной (в 2D) матрицей:

${\displaystyle u(x_0,y_0)\approx {\hat{u}}_{j,l}}$, ${\displaystyle \hat{u}=\sum _i{\hat{g}}_i⊛{\hat{k}}_i=\sum _i{\hat{A}}_i{\hat{g}}_i}$, (5)

где ${\displaystyle j,l}$ - индексы по осям ${\displaystyle oX}$ и ${\displaystyle oY}$ координат сетки ${\displaystyle \mathrm{\Omega }}$, а ${\displaystyle \hat{k}}$ - ядро свёртки, аппроксимирующее значения поля для точечного источника. Элементы ${\displaystyle \hat{k}}$ выражаются через значения ${\displaystyle k_v}$, но конкретный вид будет уточнён ниже.

По теореме о свёртке оператор ${\displaystyle {\hat{A}}_i}$ диагонализируется двумерным дискретным преобразованием Фурье. Через ${\displaystyle {ℱ}={{ℱ}}_{x,y\to v_x,v_y}}$ обозначим оператор двумерного преобразования Фурье по латерали (по осям ${\displaystyle Ox}$ и ${\displaystyle Oy}$). Тогда диагонализация ${\displaystyle {\hat{A}}_i}$ выглядит следующим образом:

${\displaystyle {\hat{A}}_i={{ℱ}}^{-1}{D}_i{ℱ}}$, (6)

где матрица ${\displaystyle D_i=Diag\left[{ℱ}{\hat{k}}_i\right]}$ - диагональная матрица собственных чисел матрицы ${\displaystyle {\hat{A}}_i}$ (прямой задачи для горизонтального слоя). Подставив (6) вместо внутренней суммы в (5), получаем численный метод решения прямой задачи:

${\displaystyle \hat{u}=\sum _{i=1..n_z}{{ℱ}}^{-1}Diag\left[{ℱ}{\hat{k}}_i\right]({ℱ}{\hat{g}}_i)}$, (7)

где ${\displaystyle \hat{u}}$ - вектор значений поля в 2D сетке, ${\displaystyle {\hat{g}}_i}$ - вектор значений физического параметра в ${\displaystyle i}$-ом горизонтальном срезе среды, а ${\displaystyle {\hat{k}}_i}$ - ${\displaystyle i}$-ое ядро дискретной свёртки. При реализации через БПФ (быстрые преобразования Фурье) вариант (7) значительно быстрее реализации (4).

Ядро свёртки

Рассмотрим свёртку в (5):

${\displaystyle {\hat{u}}_{i,j_0,l_0}={({\hat{g}}_i⊛{\hat{k}}_i)}_{j_0,l_0}={\left(\sum _{j,l}{\hat{g}}_{i,j,l}{\hat{k}}_{i,(j-j_0){modn}_x,(l-l_0){modn}_y}\right)}_{j_0,l_0}}$. (8)

Ядро ${\displaystyle {\hat{k}}_i}$ периодической дискретной свёртки (аппроксимации) не является строго равным ядру непериодической свёртки (аппроксиманту). Оно выбирается так, чтобы для центральной ячейки свёртка периодическая и непериодическая были равны. То есть, ядро периодической свёртки со смещением ${\displaystyle j_0=\lfloor n_x/2\rfloor ,l_0=\lfloor n_y/2\rfloor }$ равняется ядру непериодической с таким же смещением. В этом случае посередине сетки ${\displaystyle \mathrm{\Omega }}$ краевые эффекты в ${\displaystyle \hat{u}}$ минимальны.

Краевые эффекты

Ограниченный объём 3D модели предполагает, что интегрирование за его границами (по латерали) не осуществляется. Это утверждение эквивалентно такому случаю, при котором модель по латерали простирается до бесконечности, но в области вне модели равна нулю. В реальности так почти никогда не бывает, и поэтому значения расчётного поля около границ не соответствуют реальности даже когда значения физического параметра внутри объёма моделирования известны с высокой точностью. Более того, дискретные преобразования Фурье предполагают, что матрицы, составляющие оператор прямой задачи ${\displaystyle A_i}$, - циркулятные (в 2D), но на самом деле это не так, и расчёт через ${\displaystyle {\hat{A}}_i}$ продуцирует дополнительные краевые эффекты.

Для обработки граничных случаев используется механизм экстраполяции. Расчёт без экстраполяции называется экстраполяцией "повтор", потому что в этом случае расчёт эквивалентен случаю непериодической экстраполяции (с матрицей Тёплица), когда данные распространяются в область экстраполяции методом тайлинга (мозаичного размешения повторяющихся копий), а ядро - это ядро непериодической свёртки ${\displaystyle {\hat{k}}_i}$, дополненное нулями. Формально, экстраполированные данные ${\displaystyle {\tilde{g}}_{i,j,l}={\hat{g}}_{i,j{modn}_x,l{modn}_y}}$, а ядро свёртки ${\displaystyle {\tilde{k}}_i={\hat{k}}_i}$ для тех ${\displaystyle i}$, где определено ${\displaystyle {\hat{k}}_i}$, и 0 во всех других.

Для подавления краевых эффектов можно использовать экстраполяцию "отражением". В качестве альтернативы можно пользоваться более сложными методами экстраполяции, реализованными в утилите Экстраполяции поля с соединением противоположных краёв. При этом, если уже была произведена экстраполяция с соединением краёв, то следует выбирать "повтор".

Отличия от классического метода Приезжева

Метод основан на известном методе решения прямых задач Приезжева, основанном на преобразованиях Фурье, с важными отличиями:

1. Обработка краевых эффектов производится отражением, что показывает чуть более качественные результаты, чем реализованная Приезжевым экстраполяция линейной интерполяцией между краями сетки; также поддерживается Экстраполяция поля.
2. Есть возможность использовать гравитационное поле от прямоугольной призмы, в то время как в реализации Приезжева используется только точка.

Поле от ячейки

В рамках ГИС ИНТЕГРО полагается, что сеточной модели ставится в соответствие следующая конечноэлементная модель: для каждой точки сеточной модели строится ячейка - прямоугольная призма с центром в этой точке, а физический параметр в рамках ячейки одинаков и равен заданному в точке значению. Для гравитационной задачи есть возможность задавать три различных способа аппроксимации:

• классический расчёт поля от прямоугольной призмы (через 16 логарифмов и 8 арктангенсов);
• простая аппроксимация полем от точечной массы;
• аппроксимация прямоугольной призмы таким способом, при котором поле вдали от ячейки рассчитывается как поле от точечной массы, иначе - классический расчёт.

В случаях, когда шаг по вертикали намного меньше шага по латерали, для ячеек, находящихся близко к полю, аппроксимация полем точечной массы может иметь неприемлемо большую ошибку, но "честный" расчёт поля прямоугольной призмы может сильно замедлять выполнение алгоритма. Рекомендуется:

1. Выбирать компромиссный вариант.
2. Поле точечной массы следует выбирать для совместимости с другими системами, которые могут использовать такую аппроксимацию.

Параметры

• Исходная ТОС и исходное свойство: ТОС и имя свойства, содержащие распределение физического параметра.
• Высота поля (Z) - высота, на которой находится плоскость, содержащая сетку, в которой рассчитывается поле, м.
• Экстраполяция - Повтор или Отражение. Повтор: значение точек за краем данных равно значению точек с противоположного края. Отражение: значение точек за краем данных равно значению точек с этого же края с отступом внутрь, равным отступу наружу. Если по каким-то причинам требуется использовать стандартную экстраполяцию, то её необходимо сделать отдельно.
• Целевая ТОС и целевое свойство: ТОС и имя свойства, в которые будет записано результирующее поле. Если ТОС уже существует и совместима по привязке с исходной ТОС, то в неё будет добавлено новое свойство. Если свойство с указанным именем уже существует, будет запрошено подтверждение за перезапись свойства. Предлагается несколько имён свойств на выбор, но можно ввести собственный вариант имени свойства.